knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
# library(sjPlot)
library(multcompView)
library(modelbased)
library(hrbrthemes)
library(patchwork)
library(here)
library(readxl)
library(lmerTest)
library(blme)
# oyster length/weight data
sizdat <- read.csv(here('data/raw', 'field_data.csv'), stringsAsFactors = F) %>%
select(Exposure, Organism, Full_ID, Station, Group, Individual, Length_mm, Final_Weight) %>%
mutate(
Station = factor(Station, levels = c('OUT', 'EDG', 'MID')),
Exposure = factor(Exposure, levels = c('Pre', 'Post'))
) %>%
rename(
Weight = Final_Weight
) %>%
filter(!Organism %in% 'Mussels') %>%
gather('var', 'val', Length_mm, Weight)
# oyster dissolution data first sheet
dissdat1 <- read.csv(here('data/raw', 'Scoring data_RM_NB (1).csv')) %>%
rename(Exposure = Station) %>%
mutate(
Exposure = factor(Exposure, levels = c('OUT', 'EDG', 'MID'), labels = c('Out', 'Edge', 'Mid'))
)
# oyster dissolution data second sheet
dissdat2 <- read_excel(here('data/raw/Kelp Forest SEM ScoringNB.xlsx')) %>%
rename(Exposure = Station) %>%
mutate(
Exposure = factor(Exposure, levels = c('OUT', 'EDG', 'MID'), labels = c('Out', 'Edge', 'Mid'))
) %>%
rename(`Average.Score` = `Avg Score`) %>%
filter(!is.na(`Average.Score`))
# combine the two files
# make a unique replicate columns for two replicates at each location/species
dissdat <- bind_rows(dissdat1, dissdat2) %>%
unite('newgrp', Exposure, Group, remove = F) %>%
mutate(newgrp = factor(newgrp, labels = 1:length(unique(newgrp)))) %>%
select(-Rep)
# size models
sizmod <- sizdat %>%
mutate(
Station = case_when(
Station %in% c('EDG', 'OUT') ~ 'Edge/Out',
T ~ 'Mid'
)
) %>%
group_by(Organism, var) %>%
nest() %>%
mutate(
lmmod = map(data, ~aov(val ~ Exposure * Station, data = .x)),
phmod = map(lmmod, TukeyHSD),
ltmod = map(phmod, function(x){
modhsd <- x %>%
.$`Exposure:Station` %>%
data.frame
pvals <- modhsd$p.adj
names(pvals) <- rownames(modhsd)
lets <- multcompLetters(pvals)
return(lets)
}),
plors = pmap(list(Organism, var, lmmod), function(Organism, var, lmmod){
sjPlot::plot_model(lmmod, type = 'int', mdrt.values = 'meansd') + ggtitle(paste(Organism, var, sep = ', '))
})
)
# dissolution models
dissmod <- dissdat %>%
group_by(Organism) %>%
nest %>%
mutate(
lmmod = pmap(list(data), function(data){
# out <- lm(Average.Score ~ Exposure, data = data)
tomod <- data %>%
mutate(
newgrp = fct_drop(newgrp),
newgrp = as.numeric(newgrp)
)
# out <- blmer(Average.Score ~ Exposure + (1|newgrp:Exposure), data = tomod)
# out <- lmer(Average.Score ~ Exposure + (1|newgrp), data = tomod, REML = F) # this is the same as above if all newgrp were unique
# out <- lmer(Average.Score ~ Exposure + (1|Group:Exposure), data = tomod, REML = F) # this is also equivalent
out <- lm(Average.Score ~ Exposure, data = tomod)
return(out)
}),
anomod = map(lmmod, function(x){
out <- anova(x)
return(out)
}),
plomod = pmap(list(Organism, lmmod), function(Organism, lmmod){
# estimates
mnsval <- estimate_means(lmmod, 'Exposure')
cnsval <- estimate_contrasts(lmmod, 'Exposure') %>%
mutate(
sig = ifelse(p < 0.05, 'sig', 'notsig'),
sig = factor(sig,levels = c('notsig', 'sig'), labels = c(' not significant', 'significant'))
) %>%
unite('Contrast', Level1, Level2, sep = '-')
# sample size
n <- lmmod$model %>% nrow
# n <- lmmod@frame %>% nrow
# labels
subttl <- paste0(Organism, ' oyster')
ttl1 <- NULL
ttl2 <- NULL
ylab1 <- 'Estimated means (+/- 95% CI)'
ylab2 <- 'Estimated differences (+/- 95% CI)'
captns <- NULL#paste0('Significance is where CI does not include zero, alpha = 0.05, total n = ', n)
if(Organism == 'Pacific'){
ttl1 <- '(a) Treatment estimates'
ttl2 <- '(b) Treatment differences'
captns <- NULL
# ylab1 <- NULL
# ylab2 <- NULL
}
# mean esimate plots
p1 <- ggplot(mnsval, aes(x = Exposure, y = Mean)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = CI_low, ymax = CI_high), colour = 'black', linewidth = 1) +
labs(x = NULL, y = ylab1, title = ttl1, subtitle = subttl) +
theme_ipsum() +
theme(
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank()
) +
coord_flip()
# contrast plots
p2 <- ggplot(cnsval, aes(x = Contrast, y = Difference, colour = sig)) +
geom_point(aes(colour = sig), size = 3) +
geom_errorbar(aes(ymin = CI_low, ymax = CI_high, colour = sig), linewidth = 1) +
labs(x = NULL, y = ylab2, title = ttl2, subtitle = subttl,
caption = captns) +
theme_ipsum() +
theme(
legend.title = element_blank(),
legend.position = 'bottom',
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank()
) +
geom_hline(yintercept = 0, linetype = 'dotted', size = 1) +
scale_colour_manual(drop = F, values = c('black', 'tomato1')) +
coord_flip()
list(p1, p2)
})
)
plofun <- function(x){
psig=as.numeric(apply(x$`Exposure:Station`[,2:3],1,prod)>=0)+1
op=par(mar=c(4.2,9,3.8,2))
plot(x,col=psig,yaxt="n")
for (j in 1:length(psig)){
axis(2,at=j,labels=rownames(x$`Exposure:Station`)[length(psig)-j+1],
las=1,cex.axis=.8,col.axis=psig[length(psig)-j+1])
}
par(op)
}
sizmod$plors[[1]]
sizmod$ltmod[[1]]
## Post:Edge/Out Pre:Mid Post:Mid Pre:Edge/Out
## "a" "b" "a" "b"
plofun(sizmod$phmod[[1]])
sizmod$plors[[2]]
sizmod$ltmod[[2]]
## Post:Edge/Out Pre:Mid Post:Mid Pre:Edge/Out
## "a" "b" "c" "b"
plofun(sizmod$phmod[[2]])
sizmod$plors[[3]]
sizmod$ltmod[[3]]
## Post:Edge/Out Pre:Mid Post:Mid Pre:Edge/Out
## "a" "b" "a" "b"
plofun(sizmod$phmod[[3]])
sizmod$plors[[4]]
sizmod$ltmod[[4]]
## Post:Edge/Out Pre:Mid Post:Mid Pre:Edge/Out
## "a" "ab" "c" "b"
plofun(sizmod$phmod[[4]])
Alternative plots, these show the post-hoc comparisons of the locations from an ANOVA of treatment x location, see the above plots for all comparisons:
# Olympia
toplo <- sizmod %>%
mutate(
phmod = purrr::map(phmod, function(x){
x$`Exposure:Station` %>%
data.frame %>%
mutate(
comp = row.names(.)
)
})
) %>%
select(Organism, var, phmod) %>%
unnest('phmod') %>%
filter(comp %in% c('Post:Edge/Out-Pre:Edge/Out', 'Post:Mid-Pre:Mid')) %>%
mutate(
comp = factor(comp,
levels = c('Post:Edge/Out-Pre:Edge/Out', 'Post:Mid-Pre:Mid'),
labels = c('Out/Edge', 'Mid')
),
pcol = ifelse(p.adj <= 0.05, 'p < 0.05', 'ns'),
pcol = factor(pcol, levels = c('p < 0.05', 'ns')),
var = factor(var, levels = c('Length_mm', 'Weight'), labels = c('Delta~Length (mm)', 'Delta~Weight (g)')),
Organism = factor(Organism, levels = c('Olympia', 'Pacific'))
)
p1 <- ggplot(toplo, aes(y = comp, x = diff, linetype = pcol, col = Organism)) +
geom_vline(xintercept = 0, color = 'black', linetype = 'dotted') +
geom_point(position = position_dodge(width = 0.3), size = 3) +
geom_errorbar(aes(xmin = lwr, xmax = upr), position = position_dodge(width = 0.3), width = 0.2) +
facet_wrap(~var, scales = 'free_x', ncol = 2, strip.position = 'bottom', labeller = label_parsed) +
scale_linetype_manual(values = c('solid', 'dashed'), drop = F) +
theme_ipsum(plot_margin = margin(5, 5, 5, 5)) +
theme(
strip.placement = 'outside',
panel.grid.minor = element_blank(),
panel.grid.major.y = element_blank(),
strip.background = element_blank(),
axis.title.y = element_blank(),
legend.title = element_blank(),
legend.position = 'top',
axis.title.x = element_text(hjust = 0.5),
strip.text.x = element_text(hjust = 0.5),
strip.text = element_text(size = 9),
panel.background = element_rect()
) +
labs(
y = 'Location',
x = 'Effect size',
subtitle = '(a) Estimated means'
)
# Olympia
toplo <- sizmod %>%
mutate(
phmod = purrr::map(phmod, function(x){
x$`Exposure:Station` %>%
data.frame %>%
mutate(
comp = row.names(.)
)
})
) %>%
select(Organism, var, phmod) %>%
unnest('phmod') %>%
mutate(
pcol = ifelse(p.adj <= 0.05, 'p < 0.05', 'ns'),
pcol = factor(pcol, levels = c('p < 0.05', 'ns')),
var = factor(var, levels = c('Length_mm', 'Weight'), labels = c('Delta~Length (mm)', 'Delta~Weight (g)')),
Organism = factor(Organism, levels = c('Olympia', 'Pacific'))
)
p2 <- ggplot(toplo, aes(y = comp, x = diff, linetype = pcol, col = Organism)) +
geom_vline(xintercept = 0, color = 'black', linetype = 'dotted') +
geom_point(position = position_dodge(width = 0.3), size = 3) +
geom_errorbar(aes(xmin = lwr, xmax = upr), position = position_dodge(width = 0.3), width = 0.2) +
facet_wrap(~var, scales = 'free_x', ncol = 2, strip.position = 'bottom', labeller = label_parsed) +
scale_linetype_manual(values = c('solid', 'dashed'), drop = F) +
theme_ipsum(plot_margin = margin(5, 5, 5, 5)) +
theme(
strip.placement = 'outside',
panel.grid.minor = element_blank(),
panel.grid.major.y = element_blank(),
strip.background = element_blank(),
axis.title.y = element_blank(),
legend.title = element_blank(),
legend.position = 'top',
axis.title.x = element_text(hjust = 0.5),
strip.text.x = element_text(hjust = 0.5),
strip.text = element_text(size = 9),
panel.background = element_rect()
) +
labs(
y = 'Location',
x = 'Effect size',
subtitle = '(b) Multiple comparisons'
)
p <- p1 + p2 + plot_layout(ncol = 1, heights = c(0.8, 1))
# pout
jpeg(here('figs/fieldcomp.jpeg'), height = 8, width = 7, family = 'serif', units = 'in', res = 500)
print(p)
dev.off()
knitr::include_graphics(here('figs/fieldcomp.jpeg'))
# size, species models
sizsppmod <- sizdat %>%
mutate(
Station = case_when(
Station %in% c('EDG', 'OUT') ~ 'Edge/Out',
T ~ 'Mid'
)
) %>%
group_by(Organism, var) %>%
mutate(val2 = as.numeric(scale(val))) %>%
pivot_longer(cols = val:val2, names_to = 'scl', values_to = 'val') %>%
group_by(var, scl) %>%
nest() %>%
mutate(
lmmod = map(data, ~aov(val ~ Exposure * Station * Organism, data = .x)),
phmod = map(lmmod, TukeyHSD),
ltmod = map(phmod, function(x){
modhsd <- x %>%
.$`Exposure:Station:Organism` %>%
data.frame
pvals <- modhsd$p.adj
names(pvals) <- rownames(modhsd)
lets <- multcompLetters(pvals)
return(lets)
})
)
toplo <- sizsppmod %>%
mutate(
phmod = purrr::map(phmod, function(x){
x$`Exposure:Station:Organism` %>%
data.frame %>%
mutate(
comp = row.names(.)
)
})
) %>%
select(var, scl, phmod) %>%
unnest('phmod') %>%
mutate(
pcol = ifelse(p.adj <= 0.05, 'p < 0.05', 'ns'),
pcol = factor(pcol, levels = c('p < 0.05', 'ns'))
)
thm <- theme_ipsum() +
theme(
strip.placement = 'outside',
panel.grid.minor = element_blank(),
panel.grid.major.y = element_blank(),
strip.background = element_blank(),
axis.title.y = element_blank(),
legend.title = element_blank(),
legend.position = 'top',
axis.title.x = element_text(hjust = 0.5),
strip.text.x = element_text(hjust = 0.5),
strip.text = element_text(size = 9),
panel.background = element_rect()
)
toplo1 <- toplo %>%
filter(scl == 'val') %>%
mutate(
var = factor(var, levels = c('Length_mm', 'Weight'), labels = c('Delta~Length (mm)', 'Delta~Weight (g)'))
)
p1 <- ggplot(toplo1, aes(y = comp, x = diff, linetype = pcol)) +
geom_vline(xintercept = 0, color = 'black', linetype = 'dotted') +
geom_point(size = 3) +
geom_errorbar(aes(xmin = lwr, xmax = upr), width = 0.2) +
facet_wrap(~var, scales = 'free_x', ncol = 2, strip.position = 'bottom', labeller = label_parsed) +
scale_linetype_manual(values = c('solid', 'dashed'), drop = F) +
thm +
labs(
y = 'Location',
x = 'Effect size'
)
toplo2 <- toplo %>%
filter(scl == 'val2') %>%
mutate(
var = factor(var, levels = c('Length_mm', 'Weight'), labels = c('Delta~Length (scaled)', 'Delta~Weight (scaled)')),
)
p2 <- ggplot(toplo2, aes(y = comp, x = diff, linetype = pcol)) +
geom_vline(xintercept = 0, color = 'black', linetype = 'dotted') +
geom_point(size = 3) +
geom_errorbar(aes(xmin = lwr, xmax = upr), width = 0.2) +
facet_wrap(~var, scales = 'free_x', ncol = 2, strip.position = 'bottom', labeller = label_parsed) +
scale_linetype_manual(values = c('solid', 'dashed'), drop = F) +
thm +
labs(
y = 'Location',
x = 'Effect size'
)
# filter scaled post-hoc comparisons to those that are relevant
toplo3 <- toplo2 %>%
filter(comp %in% c("Pre:Mid:Pacific-Pre:Mid:Olympia", "Pre:Mid:Pacific-Post:Edge/Out:Pacific",
"Post:Mid:Pacific-Pre:Mid:Pacific", "Post:Mid:Pacific-Pre:Edge/Out:Pacific",
"Post:Mid:Pacific-Post:Mid:Olympia", "Post:Mid:Pacific-Post:Edge/Out:Olympia",
"Post:Edge/Out:Pacific-Post:Mid:Olympia", "Post:Edge/Out:Pacific-Post:Edge/Out:Olympia")
)
p3 <- ggplot(toplo3, aes(y = comp, x = diff, linetype = pcol)) +
geom_vline(xintercept = 0, color = 'black', linetype = 'dotted') +
geom_point(size = 3) +
geom_errorbar(aes(xmin = lwr, xmax = upr), width = 0.2) +
facet_wrap(~var, scales = 'free_x', ncol = 2, strip.position = 'bottom', labeller = label_parsed) +
scale_linetype_manual(values = c('solid', 'dashed'), drop = F) +
thm +
labs(
y = 'Location',
x = 'Effect size'
)
# pout
jpeg(here('figs/fieldcomp2.jpeg'), height = 7, width = 8.5, family = 'serif', units = 'in', res = 500)
print(p1)
dev.off()
# pout
jpeg(here('figs/fieldcomp3.jpeg'), height = 7, width = 8.5, family = 'serif', units = 'in', res = 500)
print(p2)
dev.off()
# pout
jpeg(here('figs/fieldcomp4.jpeg'), height = 5, width = 8.5, family = 'serif', units = 'in', res = 500)
print(p3)
dev.off()
knitr::include_graphics(here('figs/fieldcomp2.jpeg'))
knitr::include_graphics(here('figs/fieldcomp3.jpeg'))
knitr::include_graphics(here('figs/fieldcomp4.jpeg'))
sclex <- sizsppmod %>%
select(-lmmod, -phmod, -ltmod) %>%
unnest('data') %>%
mutate(
var1 = factor(var, levels = c('Length_mm', 'Weight'), labels = c('Length (mm)', 'Weight (g)')),
var2 = factor(var, levels = c('Length_mm', 'Weight'), labels = c('Length (scaled)', 'Weight (scaled)'))
)
thm <- theme_ipsum() +
theme(
strip.placement = 'outside',
panel.grid.minor = element_blank(),
panel.grid.major.y = element_blank(),
strip.background = element_blank(),
legend.title = element_blank(),
legend.position = 'top',
axis.title.x = element_text(hjust = 0.5),
axis.title.y = element_text(hjust = 0.5, size = 12),
strip.text.x = element_text(hjust = 0.5),
strip.text.y = element_text(hjust = 0.5),
panel.spacing = unit(0.5, 'lines')
)
toplo1 <- sclex %>%
filter(scl == 'val')
p1 <- ggplot(toplo1, aes(x = val, fill = Organism)) +
geom_density(alpha = 0.7) +
facet_grid(Exposure ~ var1, switch = 'x', scales = 'free_x') +
thm +
labs(
x = NULL,
y = 'Density',
subtitle = 'Measured values'
)
toplo2 <- sclex %>%
filter(scl == 'val2')
p2 <- ggplot(toplo2, aes(x = val, fill = Organism)) +
geom_density(alpha = 0.7) +
facet_grid(Exposure ~ var2, switch = 'x', scales = 'free_x') +
thm +
labs(
x = NULL,
y = 'Density',
subtitle = 'Scaled values'
)
# pout
jpeg(here('figs/fieldcomp5.jpeg'), height = 7, width = 8.5, family = 'serif', units = 'in', res = 500)
print(p1)
dev.off()
jpeg(here('figs/fieldcomp6.jpeg'), height = 7, width = 8.5, family = 'serif', units = 'in', res = 500)
print(p2)
dev.off()
knitr::include_graphics(here('figs/fieldcomp5.jpeg'))
knitr::include_graphics(here('figs/fieldcomp6.jpeg'))
p1 <- dissmod$plomod[[1]]
p2 <- dissmod$plomod[[2]]
p <- p2[[1]] + p2[[2]] + p1[[1]] + p1[[2]] + plot_layout(ncol = 2, guides = 'collect') &
theme(
legend.position = 'bottom',
plot.margin = unit(c(0.1, 0.1, 0.1, 0.1), 'in')
)
jpeg(here('figs/fielddiss.jpeg'), height = 7, width = 7, family = 'serif', units = 'in', res = 500)
print(p)
dev.off()
knitr::include_graphics(here('figs/fielddiss.jpeg'))
Pacific ANOVA:
dissmod %>% filter(Organism == 'Pacific') %>% pull(anomod) %>% .[[1]]
## Analysis of Variance Table
##
## Response: Average.Score
## Df Sum Sq Mean Sq F value Pr(>F)
## Exposure 2 13.481 6.7406 13.604 3.295e-05 ***
## Residuals 39 19.323 0.4955
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Pacific contrasts:
dissmod %>% filter(Organism == 'Pacific') %>% pull(lmmod) %>% .[[1]] %>% estimate_contrasts('Exposure')
## Marginal Contrasts Analysis
##
## Level1 | Level2 | Difference | 95% CI | SE | t(39) | p
## --------------------------------------------------------------------
## Edge | Mid | 0.99 | [ 0.29, 1.69] | 0.28 | 3.54 | 0.002
## Out | Edge | 0.24 | [-0.48, 0.96] | 0.29 | 0.85 | 0.402
## Out | Mid | 1.24 | [ 0.61, 1.86] | 0.25 | 4.96 | < .001
##
## Marginal contrasts estimated at Exposure
## p-value adjustment method: Holm (1979)
Olympia ANOVA:
dissmod %>% filter(Organism == 'Olympia') %>% pull(anomod) %>% .[[1]]
## Analysis of Variance Table
##
## Response: Average.Score
## Df Sum Sq Mean Sq F value Pr(>F)
## Exposure 2 10.439 5.2194 11.362 0.0001423 ***
## Residuals 37 16.997 0.4594
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Olympia contrasts:
dissmod %>% filter(Organism == 'Olympia') %>% pull(lmmod) %>% .[[1]] %>% estimate_contrasts('Exposure')
## Marginal Contrasts Analysis
##
## Level1 | Level2 | Difference | 95% CI | SE | t(37) | p
## --------------------------------------------------------------------
## Edge | Mid | 1.25 | [ 0.56, 1.95] | 0.28 | 4.52 | < .001
## Out | Edge | -0.42 | [-1.11, 0.28] | 0.28 | -1.51 | 0.141
## Out | Mid | 0.83 | [ 0.21, 1.46] | 0.25 | 3.37 | 0.004
##
## Marginal contrasts estimated at Exposure
## p-value adjustment method: Holm (1979)